The dynamics of developing Cooper pairing at finite temperatures 
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We study the time evolution of a system of fermions with pairing interactions at a finite tempera- 
ture. The dynamics is triggered by an abrupt increase of the BCS coupling constant. We show that 
if initially the fermions are in a normal phase, the amplitude of the BCS order parameter averaged 
over the Boltzman distribution of initial states exhibits damped oscillations with a relatively short 
decay time. The latter is determined by the temperature, the single-particle level spacing, and the 
ground state value of the BCS gap for the new coupling. In contrast, the decay is essentially absent 
when the system was in a superfluid phase before the coupling increase. 



Considerable progress has been made over the past few 
years in understanding the dynamical fermionic pairing 
in response to fast perturbations [lRa . Recent interest 
in this long-standing problem [10l4l2| has been motivated 
by experiments on cold atomic fermions with tunable 
interactions IH 13 . even though other systems have also 
been considered 15L 

El- 

The general picture that emerged from the theory work 
is that as a result of the perturbation, e.g. a sudden 
change of the coupling constant, the system of fermions 
with pairing interactions can reach a variety of dynamical 
phases with properties quite distinct from the equilibrium 
ones [7H9(. For example, a steady state characterized 
by undamped periodic oscillations of the time dependent 
Bardeen-Cooper-Schrieffer (BCS) order parameter A(£) 
and a gapless steady stated, 0], A(£) = 0, have been 
identified. 

Periodic oscillations occur in particular when at £ = 
the fermions are described by a many-body wave function 
with a seed gap A^ much smaller than the ground state 
gap Ao- As a result of the Cooper instability of the initial 
state the order parameter starts to grow exponentially, 
A(£) = Aie Aot , and reaches the ground state value at 
time r/2 = ln(Ao/Ai)/Ao- However, in the absence of 
the energy relaxation the system does not equilibrate and 
it can be shown that |A(£)| is periodic in time with a 
period r[l|,[i|. 

In this Letter we study the effect of temperature fluc- 
tuations on the non-adiabatic dynamics of fermions with 
attractive interaction Q. Suppose the system is initially 
in equilibrium at a finite temperature T. At £ = 
the dynamics is triggered by an abrupt increase of the 
pairing strength and a certain quantity is measured at 
a later time. This process is repeated many times for 
each data point as is typical for measurements in atomic 



rameter, (|A(t)|), displays exponentially damped oscilla- 
tions with a decay time (see also Fig. [IJ 



gases [171 UM- ^e are therefore interested in dynamical 



quantities averaged over the Boltzman distribution of ini- 
tial states. 

Our main results are as follows. We show that, if be- 
fore the coupling increase the system is in a normal phase 
at temperature T, the average amplitude of the order pa- 
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where (r) is the average oscillation period and 8 is the 
single particle level spacing. Here and below we assume 
5<T< Ao. Expression {1} is accurate up to a prefac- 
tor of order one under the logarithm. 

For typical values for cold atomic fermions 15J [14 1 
Eq. |T]) yields io/( r ) ~ 1 — 3, i.e. there are only a few 
regular oscillations before the dephasing sets in. In con- 
trast, for the paired initial phase, we demonstrate that 
to oc 1/ V5 indicating that the decay time effectively di- 
verges as the temperature is decreased below the critical 
temperature of the initial phase. 

We emphasize that each time the coupling is switched 
a particular initial condition is selected and the sys- 
tem goes into a state with periodic |A(£)|. However, 
whether the oscillations are seen in an ensemble averaged 
measurement depends on the quantity being measured. 
For example, it seems difficult to observe many of them 
in (|A(£)|). On the other hand, since the fluctuations 
of the oscillation frequency are smallQ (see also below 
Eq. (|13p ). it can in principle be obtained e.g. from the en- 
semble averaged radio frequency absorption spectra[l9(. 

The decay time Q]) can be qualitatively understood 
as follows. In the normal state a nonzero initial value 
of the order parameter Aj is due to fluctuations, which 
in mesoscopic samples are governed by an energy scale 



y/TS [2CJ, |21|. Changing Aj by a factor of order one 
in the expression for the period Aor = 21n(Ao/Aj) 
leads to changes in the period AqSt ~ 1. Then, one 
expects the average of |A(£)| over all possible values 
of A, to dephase after t/St oscillations, i.e. on to oc 
ln 2 (Ao/v / T^)/Ao timescale. Note that the average pe- 
riod (r) ex ln(Ao/v / 7^)/A and oscillation frequency re- 
main finite. In the superfluid state the order parame- 
ter has a macroscopic thermal average Aj, while typical 
thermal fluctuations VT6 -C A-j. In this case repeating 
the above argument, we obtain AqSt ~ ^/TS/Ai and 
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Ao^o ~ A, ln 2 (Ao/Ai)/\/rJ, i.e. an extremely long de- 
cay time. 

The non-stationary Cooper pairing at times much 
shorter than the energy relaxation time can be described 
by the BCS model 



H 



XS^ct^c^CkiCkf, (2) 
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where 6j are the single fermion energies relative to the 
Fermi level, S is the mean spacing between ey, A is the 
dimensionless BCS coupling constant, and dj a are the 
fermionic annihilation operators. 

In the time-dependent BCS mean-field approach[l| the 
many-body wave function is a product state 

!*(*))= LI ( U m(t)+V m (t)cl/ ml )\0), (3) 
n m =0,2 

where u rn {t) and v m (t) are the Bogoliubov amplitudes 
and the product is taken only over unoccupied (n m = 0) 
and doubly occupied (n m — 2) levels. Singly occupied 
levels are excluded since their occupation numbers are 
conserved by the Hamiltonian ([2]). 

The time evolution of the system is governed by the 
Bogoliubov-de Gennes equations 
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where A = A<5 



^ m ^rn^m- These equations can be cast 
into the form of equations of motion for classical spins [l| 

2b m x s m , b m — ( Aje, Ay, c m ) , (5) 

where A x and — A y are the real and imaginary parts of 
A = XS J2 m s~ and the components of spins are related 
to Bogoliubov amplitudes u m and v m as follows 
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For example, according to Eqs. (|3I6[) the Fermi ground 
state where all states below the Fermi level are occu- 
pied and states above are empty corresponds to = 
-sgn e m /2 and s~ = 0. 

Remarkably, nonlinear systems (0| and |(SJ| turn out to 
be mtegrable[4(. The solution for A(f) can be obtained 
with the help of the Lax vector technique by intro- 
ducing 



L(„, . - 
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where w is an auxiliary parameter and z is a unit vec- 
tor along the z axis. The square of the Lax vector is 
conserved by Eq. ([5]) for any w and therefore the roots 
of L 2 (u>) = are integrals of motion. Further, one can 
show that the majority of the roots lie on continuous 
lines, while the remaining isolated roots uniquely deter- 
mine the form of |A(t)| at times t 3> 1/Aopj. For in- 
stance, for initial states close to the Fermi ground state 



there are two isolated roots, w\ = 171 and u>2 = 172, 
in the upper half plane of complex w. In this case the 
solution of Eq. © is known to be [313111311 



|A(f)| =A + dn(A + (t-r/2),fe), k 2 = 1 - 
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where dn is the Jacobi elliptic function with modulus fc, 
t is its period, and A± = I71 ± 72 1. Eq. {SJ describes 
periodic in time |A(t)| whose period and amplitude are 
controlled by A±. 

First, consider a Fermi gas at a temperature T and 
zero BCS coupling constant. At t = the coupling is 
suddenly turned on so that A ^> T, where A is the 
ground state gap for the new coupling. Before the in- 
teraction switch on the system can be in any eigenstate 
of the free Fermi gas with the probability given by its 
Boltzman weight. These states thus provide an ensemble 
of initial conditions for equations of motion ([5]) and our 
task is to evaluate the average of |A(i)| over all possible 
initial states. 

In the non-interacting problem amplitudes (u m ,v m ) 
take values (0,0), (1,0), and (0,1) corresponding to oc- 
cupancies n m — 1, 0, and 2, respectively. Note that 
they are always correlated so that s~ = u m Vm = and 
sf n = ±1/2 or indicating that the eigenstates of the 
free Fermi gas are (unstable) stationary states for the 
mean- field equations of motion (|4l5p . However, for any 
nonzero coupling they are not exact stationary states of 
the quantum Hamiltonian {2j before the mean-field de- 
coupling of the interaction term. These quantum effects 
facilitate the development of the Cooper instability and 
after a short time states of the form ([3]) with finite u m v^ 
can be used. In the spin language, the spins s TO acquire 
nonzero s~ , i.e. nonzero components in the xy plane. 

As argued in Ref. |y only spins at energies \e m \ < T <C 
Ao initially have appreciable xy components (see below) . 
It follows from Eq. ((7j) that L 2 (u>) has two isolated roots 
in the upper half plane of complex w and the order pa- 
rameter is described by Eq. l[8|). where the parameters 
A± are 
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The values of s~ are random with a distribution de- 
termined by the Boltzman distribution of initial states 
and the quantum effects discussed above. On the other 
hand, there is a large number N ~ T/S of random com- 
plex numbers in the sum ^ and as noted in Ref. 
(see Eq. (46) therein) one therefore expects the Rayleigh 
distribution 1221 



p(A_) = CA_ exp 



aA 2 _ 
4T5 
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independent of the details of the distribution of s~ . Here 
^/TS is a characteristic scale of fluctuations of A_, a is 
of order one, and C is a normalization constant. 
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Thus, averaging |A(f)| over Boltzman distributed ini- 
tial states reduces to integrating Eq. JSj) with respect to 
A_ with distribution ([TT)1), i.e. 



0.7 



Ao 



/ dn(A (t-r/2),fc)p(A_)dA_ 
Jo 
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Note that the Jacobi function dn depends on A_ through 
its modulus k = 1 — A^/Aq. For example, its period for 
A_ < A is M 



A 
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Using Eqs. (|12ll0p . we evaluate the average oscillation 
period and its standard deviation (see also Ref. @), 
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up to a factor of order one under the logarithm. The 
average frequency and its deviation are (u>) = 2~k/(t) 
and Su>/(lo) — St j (t). 

The asymptotic behavior of integral (fTTj) at large times 
t to can be evaluated using the saddle point method 
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where t is given by Eq. ([T]). We see that on the t 
time scale (|A(f)|) exponentially approaches a constant 
value smaller than the ground state gap A by a large 
factor 1ii(q;Aq/T5) Jit. The approach is oscillatory with a 
typical period close to the ensemble averaged period (r) . 

Next, we present several alternative systematic deriva- 
tions of Eq. (fT0|) and show that it is independent of 
the details of initial state distribution. First, note that 
Eqs. ([5]) are equations of motion for classical spin Hamil- 
tonian H = J2 m 2e m s* n - A££ m ,„ Sm 3 ™- As discussed 
above, before the interaction switch on the spins are 
along the z axis. Their z components take values = 
±1/2 or with independent probabilities proportional 
to the corresponding Boltzman weig ht e - 2e - s ™/ T . This 
presents a technical difficulty, since these spin configura- 
tions are (unstable) equilibria for Eqs. ([5]). 

One way to circumvent this problem is to replace 
the above ensemble of initial spin configurations with 
the Boltzman distribution of classical spins of length 
s m = 1/2. Then, each spin s m is characterized by po- 
lar and azimuthal angles 9 m and </> m with independent 
probability proportional to exp(— e m cos 6 m /T), i.e. spins 
at \e m \ < T acquire finite components in the xy plane. 
Using this probability distribution and Eq. ©, we eval- 
uate p(A_). The calculation results in Eq. (fTU|) with 
a = 2/ln(A /T) and we obtain Eqs. (fTTTB"]) and CO 




FIG. 1: Time evolution of the amplitude of the BCS order 
parameter |A(t)| averaged over initial states at temperature 
T. Numerical simulation of Eq. (Jsj) for 10 4 spins averaged 
over 10 4 realizations of initial conditions (|15|l (solid curve) is 
compared to expression (dotted curve). The time is in 
units of the decay time to, Eq. CfJ; the ground state gap is 
A = 2 x 10 3 c5 and T = 4005, where S is the level spacing. 



A distribution of the form (TO]) for \A(t = 0)| was 
obtained in Ref. [a for an ensemble of initial condi- 
tions suggested in the same reference. Note that ac- 
cording to Eq. Q \u m \ 2 and \v m \ 2 represent probabil- 
ities of zero and double occupancy, respectively, of the 
level e m . In the free Fermi gas before the interaction is 
turned on their thermal averages are (|u r „ l 2 ) = n 2 ^ and 



= (1 — ?i m ) 2 , where 



( e e m /T + i s the 



Fermi function. Averaging Eq. ((4]) with respect to e m in 
a narrow window of energies, we can replace u m and i?„ 



with 



and 



'(l-n m ) 2 , 



where 4> m is a random relative phase. Since the total 
energy of the free gas does not depend on m , they are 
assumed to have independent uniform distributions. Fur- 
ther, assuming (umV^) w (u m ) (v!^) , and using Ec 
we obtain the following initial spin configurations 



1 



= — tanh , , 
2 V2TV 
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Using Eqs. (fl~5l 15)) and uniform distributions for <fi m , we 
derive Eqs. fTTTB")) and Eq. CTO|) with a = 6, see also Fig.Q] 
Finally, Eqs. (|H13p and (fl"4"|) can be derived starting 
from the Ginzburg-Landau free energy. The advantage 
of this approach is that we can consider initial states 
with nonzero BCS coupling that is suddenly increased at 
t = 0. The ground state gap for the new coupling, Ao, is 
assumed to be much larger than that for the old coupling. 
Then, the equation L 2 (w) — has two isolated roots with 
Imu; > and the evolution of the order parameter is 
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described by Eq. (JSJ) as before. With the help of Eq. 0, 
we obtain A + w Ao and A_ w 2Aj ln(Ao/Aj), where A, 
is the gap for the old coupling. 

First, consider the case T > T c , where T c is the critical 
temperature for the old coupling. To calculate the aver- 
age of |A(t)| over initial states, we need the probability 
distribution of A_ or equivalcntly the distribution of pos- 
sible values of the gap A, before the coupling increase. 
We assume the latter is of the Ginzburg-Landau form 
Aj exp(-F(Aj)/T), where the free energy for T > T c is 



F (Aj) = ln(T/r c )|Ai| 2 /5 [20|,l2l|. Using this distribu- 
tion function and the above expressions for A± in terms 
of Aj, we again obtain Eq. I|14p. where now 
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This expression holds for T - T c > \fTJ>. 

Below the critical temperature, for T c > T c — T ^> 
y/T c S, we keep the quartic term in F(Ai) and expand 
Eq. in Aj — Aj, where Aj is the thermal average of 
the order parameter before the coupling change. Using 
a saddle point method, we obtain a Gaussian decay to a 
constant value on t cx A 2 ln 2 (A /A i )/\/r 5 5 timescale. 
On the other hand, the dynamics at times this long is 
likely not described by the Hamiltonian © that docs 
not account for energy relaxation. Thus, we see that 
the dephasing of ensemble averaged oscillations due to 
thermal fluctuations is effectively absent when the dy- 
namics is started in the paired phase. The reason is that 
in this case the order parameter has a macroscopic ini- 
tial average much larger than its thermal fluctuations. 
The fast dephasing above T c crosses over into a slow de- 
phasing below T c in a narrow window of temperatures 
\T-T C \~ VTJ. 

In conclusion, we studied the effect of thermal fluctu- 
ations on the dynamics of fermions with pairing inter- 
actions triggered by an abrupt increase of the pairing 
strength. We showed that if the system is in the nor- 
mal phase before the coupling increase, the amplitude of 
the order parameter averaged over the Boltzman distri- 
bution of initial states exhibits damped oscillations with 
relatively short decay time ([Tjl. see Eq. (|14[). On the 
other hand, the damping is essentially absent when the 
dynamics starts from the superfluid phase. 

An interesting problem is to determine the time evolu- 
tion described by the quantum Hamiltonian |2]) at T = 
starting from the Fermi ground state, i.e. the ground 
state of the Hamiltonian ((2]) for A = 0. Extending 
the above considerations to this case, one might expect 
damped oscillations due to quantum fluctuations with- 
out ensemble averaging. If this is the case, an estimate 



for the decay time can be obtained by replacing the 
temperature T in Eq. ([1]) with the level spacing 8, i.e. 
t ~ ln 2 (A /S)/A . 
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